function zSpread = bndspread(SpotRateData, Price, CouponRate, Settle, ...
    Maturity,varargin)
% BNDSPREAD : Static spread over spot curve
%
% Spread = bndspread(SpotInfo, Price, Coupon, Settle, Maturity)
%
% Spread = bndspread(SpotInfo, Price, Coupon, Settle, Maturity, ...
%       Period, Basis, EndMonthRule, IssueDate, FirstCouponDate, ...
%       LastCouponDate)
%
% Spread = bndspread(SpotInfo, Price, Coupon, Settle, Maturity,
%      'Param1','Value1',...)
%
%   Optional Inputs: Period, Basis, EndMonthRule, IssueDate, FirstCouponDate,
%                    LastCouponDate, StartDate, Face, LastCouponInterest,
%                    CompoundingFrequency, DiscountBasis
%
%     Note that optional inputs can be specified as parameter value pairs.  If
%     LastCouponInterest, CompoundingFrequency or DiscountBasis are input,
%     optional inputs must be specified as parameter value pairs.
%     Otherwise, optional inputs may be specified by order according to the
%     help.
%
% Required Inputs:
%
% SpotInfo  - [Matrix of two columns]
%             First column is the SpotDate, and second column is the
%             zero-rate corresponding to maturities on the SpotDate.
%
%             (It is highly recommended that the spot-rates are spaced as
%             evenly apart as possible, perhaps one that is built from
%             3-months deposit rates.)
%
%             Example:
%              SpotInfo = ...
%                 [datenum('2-Jan-2004') ,  0.03840;
%                  datenum('2-Jan-2005') ,  0.04512;
%                  datenum('2-Jan-2006') ,  0.05086];
%
%          Price  - [SCALAR or VECTOR]
%                   Price for every $100 notional of bonds whose spreads
%                   are to be computed
%                   Example: 105.484
%
%          Coupon - [SCALAR or VECTOR]
%                   Annual Coupon rate of bonds whose spreads are to be
%                   computed
%                   Example: 0.04375
%
%          Settle - [SCALAR or VECTOR of IDENTICAL ELEMENTS]
%                   Settlement date of bonds whose spread are to be
%                   computed.
%                   Example: datenum('26-Nov-2002')
%
%        Maturity - [SCALAR or VECTOR]
%                   Maturity date of bonds whose spread are to be computed
%                   Example: datenum('15-Oct-2006')
%
% Optional Inputs:
%
%          Period - [SCALAR or VECTOR]
%                   Number of coupon payment in a year.
%                   Default is 2 (semiannual)
%                   Possible values are 1, 2, 3, 4, 6, 12.
%
%           Basis - [SCALAR or VECTOR]
%                   Basis of bonds whose spread are to be computed.
%
%                   0 - actual/actual(default)
%                   1 - 30/360 SIA
%                   2 - actual/360
%                   3 - actual/365
%                   4 - 30/360 PSA
%                   5 - 30/360 ISDA
%                   6 - 30E/360
%                   7 - actual/365 Japanese
%                   8 - actual/actual ISMA
%                   9 - actual/360 ISMA
%                  10 - actual/365 ISMA
%                  11 - 30/360 ISMA
%                  12 - actual/365 ISDA
%                  13 - bus/252
%
%    EndMonthRule - [SCALAR or VECTOR]
%                   Use 1 for on (default) or 0 for off. Will make sure
%                   all payments are in the end of months when one of
%                   coupon is on the last day of a 30-day month.
%
%       IssueDate - [SCALAR or VECTOR]
%                   Also known as Dated date.
%                   Example: datenum('25-Oct-2001')
%
% FirstCouponDate - [SCALAR or VECTOR]
%                   Date when a bond makes its first coupon payment.
%                   When FirstCouponDate and LastCouponDate are both
%                   specified, FirstCouponDate takes precedence in
%                   determining the coupon payment structure.
%                   Example: datenum('15-April-2002')
%
%  LastCouponDate - [SCALAR or VECTOR]
%                   Last coupon date of a bond prior to the
%                   maturity date.
%
%       StartDate - Forward starting date of payments.
%
%  LastCouponInterest - Compounding convention for computing the yield of
%                        a bond in the last coupon period, i.e.: with only
%                        the last coupon and the face value to be repaid.
%                        Choices are 'simple' or 'compound'.
%
%  CompoundingFrequency - Compounding frequency for yield calculation.  By
%                          default, SIA bases (0-7) and BUS/252 use a semi-annual
%                          compounding convention and ISMA bases (8-12) use
%                          an annual compounding convention.
%
%  DiscountBasis - Basis used to compute the discount factors for
%                   computing the yield.  The default behavior is for SIA
%                   bases to use the actual/actual day count to compute
%                   discount factors, and for ISMA day counts and BUS/252
%                   to use the specified basis.
%
% Outputs:
%
%          Spread - [SCALAR or VECTOR]
%                   Static Spread to benchmark, in basis points.
%
% Example:
% % Computing a FNMA 4 3/8 spread over treasury spot-curve
%
% % Building Spot Curve
% RefMaturity = [datenum('02/27/2003');
%                datenum('05/29/2003');
%                datenum('10/31/2004');
%                datenum('11/15/2007');
%                datenum('11/15/2012');
%                datenum('02/15/2031')];
%
% RefCpn = [0;
%           0;
%           2.125;
%           3;
%           4;
%           5.375] / 100;
%
% RefPrices =  [99.6964;
%               99.3572;
%              100.3662;
%               99.4511;
%               99.4299;
%              106.5756];
%
% RefBonds = [RefPrices, RefMaturity, RefCpn];
% Settle   = datenum('26-Nov-2002');
% [ZeroRates, CurveDates] = zbtprice(RefBonds(:,2:end), RefPrices, Settle)
%
% % FNMA 4 3/8 maturing 10/06 at 4.30 pm Tuesday, Nov 26, 2002
% Price    = 105.484;
% Coupon   = 0.04375;
% Maturity = datenum('15-Oct-2006');
%
% Period = 2;
% Basis  = 1;
% SpotInfo = [CurveDates, ZeroRates];
%
% SpreadInBP = bndspread(SpotInfo, Price, Coupon, Settle, Maturity, Period,
% Basis)

% Copyright 2003-2009 The MathWorks, Inc.
% $Revision: 1.1.6.15 $  $Date: 2009/11/05 16:59:01 $

% parse inputs:
if nargin < 5
    error('Finance:bndspread:tooFewInputs',...
        'You must enter SpotInfo, Price, Coupon, Settle, and Maturity.');
end

% Check Settle to make sure its scalar
Settle = datenum(Settle);
if numel(Settle) > 1
    Settle = Settle(1);
    warning('Finance:bndspread:vectorSettle',...
        'Settle must be a scalar - first settle date has been selected')
end

% Check to see whether we have the case of ordered inputs or PV pairs
if ~isempty(varargin)
    if ischar(varargin{1})
        
        p = inputParser;
        
        p.addParamValue('discountbasis',NaN);
        p.addParamValue('compoundingfrequency',NaN);
        
        p.addParamValue('period',[]);
        p.addParamValue('basis',[]);
        p.addParamValue('endmonthrule',[]);
        p.addParamValue('issuedate', []);
        p.addParamValue('firstcoupondate', []);
        p.addParamValue('lastcoupondate', []);
        p.addParamValue('startdate', []);
        p.addParamValue('face',[]);
        
        try
            p.parse(varargin{:});
        catch ME
            newME = MException('finance:bndspread:optionalInputError',...
                'Error in input arguments');
            newME = addCause(newME,ME);
            throw(newME)
        end
        
        CompFreq = p.Results.compoundingfrequency;
        Basis =  p.Results.basis;
        Period =  p.Results.period;
        DiscountBasis = p.Results.discountbasis;
        
        try
            [CouponRate, ~, Maturity,Period,Basis] = instargbond(CouponRate, Settle, Maturity,Period, Basis);
        catch ME
            throwAsCaller(ME)
        end
    else
        CompFreq = NaN;
        DiscountBasis = NaN;
        
        try
            [CouponRate, ~, Maturity,Period,Basis] = instargbond(CouponRate, Settle, Maturity,varargin{:});
        catch ME
            throwAsCaller(ME)
        end
    end
else
    
    CompFreq = NaN;
    DiscountBasis = NaN;
    
    try
        [CouponRate, ~, Maturity,Period,Basis] = instargbond(CouponRate, Settle, Maturity);
    catch ME
        throwAsCaller(ME)
    end
end

if size(SpotRateData,2)<2
    error('Finance:bndspread:tooFewColumns',...
        'Need at least 2 columns: SpotDates and SpotRates.')
else
    CurveDates = datenum(SpotRateData(:,1));
    ZeroRates = SpotRateData(:,2);
end

% admit only sorted and unique set of forward rates
[CurveDates, m] = unique(CurveDates);
ZeroRates = ZeroRates(m);

% check if spotrate data is obsolete
if min(CurveDates) < Settle
    error('Finance:bndspread:invalidSettle',...
        'First spot rate data must be on, or after Settle.')
end

% appending the spot rate at Settle if not yet done
% it is going to be equal to the shortest maturity
% as measured from Settle.
if min(CurveDates) > Settle
    CurveDates(2:end+1) = CurveDates;
    CurveDates(1) = Settle;
    ZeroRates(2:end+1) = ZeroRates;
end

% Scalar expansion
try
    [CouponRate, Settle, Maturity, Period, Basis, Price] = ...
        finargsz('all', CouponRate, Settle, Maturity, Period, Basis,Price);
catch ME
    throwAsCaller(ME)
end

[CFlowAmounts, CFlowDates, CFTimes] = cfamounts(CouponRate, Settle, Maturity, varargin{:});

% Get spot rates at cash flow dates
CFSpotRates = interp1(CurveDates, ZeroRates, CFlowDates);

% Subtract price so Spread will make NPV = 0
CFlowAmounts(:,1) = CFlowAmounts(:,1) - Price;

% Spread is initially guessed to be close to the coupon at each period
x0 = CouponRate ./  Period;

% setting an option to suppress optimization results
fZeroOptions = optimset('fzero');
fZeroOptions = optimset(fZeroOptions,'Display','off');
fZeroOptions = optimset(fZeroOptions,'TolX',1e-12);

% Also, set options for FSOLVE -- in case FZERO does not converge
fSolveOptions = optimset('fsolve');
fSolveOptions = optimset(fSolveOptions, 'TolX', 1e-14);
fSolveOptions = optimset(fSolveOptions, 'TolFun', 1e-12);
fSolveOptions = optimset(fSolveOptions, 'Display', 'off');
fSolveOptions = optimset(fSolveOptions, 'LargeScale', 'off');

% Loop through each bond to find the spread
nBonds = length(Price);
X = NaN*ones(nBonds,1);
EXITFLAG = zeros(nBonds,1);

% Figure out the compounding frequency here
if isnan(CompFreq)
    CompFreq = 2*ones(nBonds,1);
    if isnan(DiscountBasis)
        i = isisma(Basis);
    else
        i = isisma(DiscountBasis);
    end
    CompFreq(i) = 1;
end

for bondidx = 1:nBonds
    
    objfun = @(x) nansum(CFlowAmounts(bondidx,:).*1 ./( 1 + ...
        (CFSpotRates(bondidx,:) + x)./CompFreq(bondidx)).^CFTimes(bondidx,:));
    
    % FZERO will iteratively solve for the value of OAS (X) that will make
    % the price of the bond equal to given market clean price
    [X(bondidx),~,EXITFLAG(bondidx)]= fzero(objfun, x0(bondidx), fZeroOptions);
    
    % if fzero was unable to find a zero, try using fsolve
    if (EXITFLAG(bondidx) < 0)
        [X(bondidx),~,EXITFLAG(bondidx)] = fsolve(objfun, x0(bondidx), fSolveOptions);
    end
    
end

if any(EXITFLAG <0)
    warning('Finance:bndspread:solutionConvergenceFailure',...
        'Could not solve for the spread for at least one of the bonds.\n');
end

% basis point is 10,000 times of Spread in decimal
zSpread = 10000*X;